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Abstract 

Recent ROSAT measurements show that the x-ray emission from iso- 
lated neutron stars is modulated at the stellar rotation period. To inter- 
pret these measurements, one needs precise calculations of the heat trans- 
fer through the thin insulating envelopes of neutron stars. We present 
nearly analytic models of the thermal structure of the envelopes of ul- 
tramagnetized neutron stars. Specifically, we examine the limit in which 
only the ground Landau level is filled. We use the models to estimate 
the amplitude of modulation expected from non-uniformities in the sur- 
face temperatures of strongly magnetized neutron stars. In addition, we 
estimate cooling rates for stars with fields B ~ 10 15 — I0 16 G which are 
relevant to models that invoke "magnetars" to account for soft 7-ray emis- 
sion from some repeating sources. 

1 Introduction 

Since the launch of the ROSAT satellite, our knowledge of isolated neutron stars 
has expanded into new realms. Before ROSAT, neutron stars were somewhat 
unique among astronomical objects. Although they had been observed over a 
range of energies from radio to ultra-high-energy gamma rays, and had been 
evoked to power a variety of astrophysical objects from pulsars to soft X-ray 
repeaters and gamma-ray bursts, one could not argue unequivocally that a single 
photon from the surface of a neutron star had ever been detected. 

For the first time, we have direct evidence for radiation from the surfaces of 
neutron stars. More than a dozen such sources have been detected by ROSAT 
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Table 1: Several pulsars with observed surface blackbody emission 



Pulsar 



References 



PSR J0437-4715 

PSR 0630+18 (Geminga) 

PSR 0656+14 



PSR 0833-45 (Vela) 



PSR 1929+10 



Becker & Triimpcr 1993 

Halpern & Holt 1992, Halpcrn & Rudcrman 1993 

Halpern & Wang 1997 

Finley, Ogelman & Kiziloglu 1992, 

Anderson et al. 1993, Greiveldinger et al. 1996, 

Possenti, Mereghetti & Colpi 1996 

Ogelman, Finley & Zuckcrman 1993 

Ogelman & Finley 1993, Greiveldinger et al. 1996 

Yancopoulos, Hamilton & Helfand 1994 



(e.g. Ogelman 1995), and more than ten have been fitted with spectra. The 
spectra divide the objects into two classes: 1) objects with hard spectra whose 
X-ray emission is best attributed to the magnetosphere, and 2) neutron stars 
whose soft flux is well-described by a blackbody spectrum. Table 1 lists several 
of those that fall into the second group with pertinent references. 

The X-ray spectra of these objects consist of a soft blackbody and hard 
power-law component. Generally, the surface temperature inferred from model- 
ing the blackbody is approximately what one would expect for a cooling neutron 
star at the characteristic age of the pulsar. If this interpretation is correct, the 
observations provide a direct probe of the structure of cooling neutron stars. 
Furthermore, the observations show that the thermal radiation is modulated 
at the rotation period of the pulsars. A strong magnetic field can modulate 
the thermal flux by causing the heat conduction in the outer layers of the star 
to be anisotropic. Both to translate these observations into constraints on the 
structure of neutron stars and to understand modulation of the radiation, we 
must have a detailed understanding of the insulating layers of the neutron-star 
crust, i.e. the envelope. 

We proceed in the spirit of Gudmundsson, Pethick & Epstein (1982) and 
concentrate our analysis on the thin region, the envelope, which insulates the 
bulk of the neutron star. The envelope is customarily defined to extend from 
zero density to p <~ 10 10 g/cm 3 , and its thickness (He) is of the order of tens of 
meters, very small compared to the radius of the star, R ~ 10 km. By limiting 
the analysis, we focus on how various physical processes affect the thermal 
structure of the envelope and the relationship between the core temperature and 
the flux emitted at the surface. An alternative point of view is to combine the 
envelope calculation with an estimate of the cooling rate due to neutrinos and 
the total heat capacity of the neutron star, yielding theoretical cooling curves 
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(e.g. Tsuruta et al. 1972, Nomoto & Tsuruta 1981, Glen & Sutherland 1980, Van 
Riper & Lamb 1981) 

Several authors have made much progress in understanding the properties 
of neutron star envelopes with and without magnetic fields. Gudmundsson, 
Pethick & Epstein (1982) numerically calculate the thermal structure for un- 
magnetizcd envelopes, and Hernquist & Applegate (1984) present analytic mod- 
els for the B = case. Tsuruta (1979), Glen & Sutherland (1980), Van Riper & 
Lamb (1981) and Urpin & Yakovlev (1980) calculate the luminosity observed at 
infinity as a function of the core temperature for several magnetic field strengths 
less than 10 14 G, including the zero-field case. Hernquist (1985) calculates the 
thermal structure of envelopes for B < 10 14 G for transport along the field, us- 
ing the electron conductivities of Hernquist (1984) which account for the quanti- 
zation of electron energies in the magnetic field in a rclativistic framework. We 
will use these conductivities in the present work; therefore, Hernquist (1985) 
provides a natural benchmark. 

Van Riper (1988) builds upon the Hernquist (1985) results by exploring 
various assumptions concerning the properties of the envelope at low densities 
and calculating profiles for many field strengths (B < 10 14 G) and core tem- 
peratures. Again, these calculations are limited to conduction along the field. 
Schaaf (1990), using the electron conductivities calculated in Schaaf (1988), 
calculates the thermal structure in two dimensions for B <j 10 11 G. Above a 
field strength of 10 12 G, the calculations are not considered reliable. Finally, 
Shibanov et al. (1995) present the temperature distribution as a function of mag- 
netic eolatitude for B = 10 12 G from a numerical solution to the two-dimensional 
thermal structure equation in a plane-parallel approximation. 

The current work complements the previous ones by extending the results 
to stronger field strengths (10 14 G < B < 10 16 G) in a semi-analytical fashion. 
We apply the approach of Hernquist & Applegate (1984) in the limit of a strongly 
magnetized envelope, and then justify and use the plane-parallel approximation 
to solve the two-dimensional structure equation. We derive separable thermal 
structure equations in the high and low temperature limits for both liquid and 
solid material. We calculate the thermal structure in terms of simple (although 
analytically intractable) integrals. 

The plane-parallel approximation has the second important advantage that 
the detailed field configuration separates from the thermal structure problem. 
Assuming that it is correct, we can synthesize the results for any field distribu- 
tion B(9, <p) as long as B is not too inhomogeneous on the scale of the envelope 
thickness (i.e. \B/VB\ > h E ). 

We find that the emission from a given surface element is a simple function 
of the location of the clement. Using this functional form, we derive light curves 
and time-dependent spectra including general relativistic effects. Although we 
closely follow the formalism of Page (1995), we calculate the two-dimensional 
thermal structure of the envelope and present results for several field strengths 
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and fluxes. We internally verify and justify the geometric simplification used to 
translate our results into observables. 



2 Preliminaries 

In extremely intense magnetic fields, the Landau energy (Ucob) of an electron 
will typically exceed its thermal energy. In these strong fields, the quantization 
of the electron energy determines the structure of the electron phase space and 
must be taken into account in calculating the thermodynamics of the electron 
gas. 

In what follows, we will use the dimensionless units 

hw B _ ?i|e| B 
P ' m e c 2 mfc 3 ~ 4.4 x 10 13 G' 1 ' 

and C - (2) 



m e c 2 5.9 x 10 9 K 

where p, is the chemical potential of the electron gas including the electron rest 
mass. Heyl (1998), Hernquist (1985), and Hernquist (1984) outline the tech- 
niques for calculating the thermodynamic properties of a magnetized electron 
gas, and we adopt their methodology for the remainder of the paper. 



3 The Low- Temperature, Strong-Field Regime 

We are specifically interested in the low temperature limit (r<(-l) and the 
regime in which only one Landau level is filled (( < ^2/3 + 1). For neutron 
stars with (3 ;> 1, this limit applies to the regions that most effectively insulate 
the isothermal core of the star. We will use the results of Yakovlev (1984) and 
Hernquist (1984) to calculate the thermal conduction in the liquid and solid 
phases. 



3.1 Degenerate Structure Equations 

If we assume that the pressure is supplied by the electrons alone, the general 
relativistic equations of thermal structure in the plane-parallel approximation 
assume the simple form 

dT = FYe_l f 1 _FS i y 1 
dp g s m u k \ g s pn J 

^ - ^-gtx-L^A (4 ) 



dz Ye' \ 9s pn 

where we have neglected the thickness of the envelope (He ~ 100 m) relative 
to the stellar radius (R). Here, m u is the atomic mass unit, Z and A are the 
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mean atomic number and mean atomic mass of the material, p is the density 
of the matter, S e is the entropy of the electron gas per unit volume, k is the 
thermal conductivity, and F and g s are the flux and the acceleration of gravity 
as measured at the surface, respectively. For completely ionized material, Y e is 
given by the product of Z/A and the ionized fraction. 

In the absence of a magnetic field, this plane-parallel approximation in- 
troduces errors of the order R s h E /R 2 ~ 0.6% where R s = 2GM/c 2 (Gud- 
mundsson, Pethick & Epstein 1982). To understand the potential errors of the 
plane-parallel treatment in the presence of a magnetic field, we compare the re- 
sults of Shibanov et al. (1995) with those of Schaaf (1990). Although Shibanov 
et al. (1995) use a one-dimensional approach, their results agree with those of 
the two-dimensional calculations by Schaaf (1990) for several surface temper- 
atures and a magnetic field of 10 12 G. In stronger magnetic fields, conduction 
perpendicular to the magnetic field is even less important, and we expect the 
one-dimensional method to be even more accurate. 

To estimate the errors in using the plane-parallel treatment in the presence of 
the magnetic field, we examine the thermal structure equation in two dimensions 
(Schaaf 1990) 



VF = -- 



-A. Q 



dr 



2 _ A dT dT 
r Kue s — + rni2 — 
or oO 



(5) 



d 



_ As dT K22 dT 



= 0, 



r sin 9 89 

where 9 is an angle along the surface of star, specifically the magnetic colatitude, 



1 ^ f , (6) 

and the the components of the thermal conduction tensor, where 1 

denotes the radial direction and 2 denotes the tangential direction. The com- 
ponents of k are found by rotating the tensor calculated by Hcrnquist (1984) 
and Yakovlev (1984) so that the ^-direction locally coincides with the radial 
direction. This gives 

«ii = K yy sin 2 ip + k zz cos 2 ip, (7) 
«22 = K yy cos 2 ip + k zz sin 2 ip, (8) 
«12 = «21 = {Kyy - k zz ) sin ip cos ip, (9) 

where ip is the angle between the local field direction and the radial direction, 
and K yy and k zz respectively are the components of the heat conduction tensor 
perpendicular and parallel to the direction of the magnetic field. 

For a uniformly magnetized neutron star ip = 9; for a dipolc field, cot?/; = 
2cot# (Greenstein & Hartke 1983), or more conveniently 

9 , 4cos 2 # ,„„, 

cos 2 ip = — . (10) 

3 cos 2 9 + 1 K ' 
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If we assume that the components of the thermal conduction matrix (n) are 
of the same order and take the maximum temperature gradient to be T c — T e s ~ 
T c radially over the thickness of the envelope, He, or tangentially over one 
radian, we obtain 

% ~ T c « e~^r^- ~ e~ A -R^- ~ 10 2 T C (11) 
o9 or h E 

where T c is the core temperature. We find that neglecting derivatives with 
respect to angle does not dramatically increase the error relative to the unmag- 
nctized plane-parallel case. 

However, this argument does not apply for 9 close to 7r/2 (i.e. where the 
magnetic field lines are parallel to the surface for a uniform or dipole field). 
Here, 

«22 = K zz > Kll = Kyy. (12) 

If we reexamine the error analysis for ip w tt/2 we find that relevant quantities 
to compare are 

K22^ and Kll r 2 e- 2A =— . (13) 
The tangential transport will exceed the radial transport if 

cos 2 ^<e 2A ^-^ (14) 

R K zz 

Comparing these values we find that if k zz ;> 10 4 K yy , the one-dimensional treat- 
ment will break down near ip = tt/2; otherwise, the plane-parallel treatment is 
adequate even at ip = tt/2. 

For regions where the magnetic field lines arc not nearly parallel to the 
surface, the plane-parallel approximation works well; consequently even for an 
arbitrary field geometry, because the envelope is thin, we ignore 8/ 89 terms in 
the structure equation compared to rdjdr terms and focus on radial heat flow. 
With these assumptions, we have (Schaaf 1990) 

k = K11 = «zz cos 2 ip + Kyy sin 2 ip. (15) 

In the low-temperature limit, we obtain the dimensionless equation 



(16) 



where X e is the electron Compton wavelength and the dimensionless flux is given 

( Y e ^-^-) = 7.83 x io-3^^, (it) 
\ m u g s c / A 56 14 
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where Z 26 = Z/26, A 56 - A/56, T cff , 6 = T cff /10 6 K, and g sM = .g s /10 14 cm/s 2 . 
T e ff is the effective blackbody temperature of the neutron star photosphere again 
as measured at the surface, which we take to be located at an optical depth of 
2/3. 

For Equation 16 to be separable for an arbitrary geometry, n zz and k vv must 
depend on r in the same fashion. For electron-ion scattering this is the case, 
so we can hope to find a simple analytic solution to the structure equation. 
Unfortunately, since the cross-section for electrons to scatter off of phonons 
depends explicitly on temperature, in the solid state the structure equation is 
separable only where the field is either purely radial or tangential. 

In the liquid state, we obtain the following structure equation, 



2 ; tt (3 2 <MC;/3) , sin 2 ^Za 2 



cos tp — 



+ 



3 Za 2 v /(2 — T 12tt 13 



Q el (C:(3)V( 2 "I 



dC \ m u g s c 

(18) 

where a is the fine-structure constant. In the solid state for electron-phonon 
scattering, we obtain 




MC;/3) 



Q ep (C;P) 



(19) 
(20) 



The functions <j> and Q are defined and calculated in Yakovlev (1984) and 
Hernquist (1984) and we take u_ 2 ~ 13 (Yakovlev & Urpin 1980, Potekhin 
& Yakovlev 1996) for a body-centered cubic lattice. 



4 The High Temperature Regime 

In the nondegenerate regime we assume that most of the heat is transported 
by photons and that free-free absorption provides the opacity. We take the 
unmagnetized thermal conductivity to be of the Kramer's form (Silant'ev & 
Yakovlev 1980), 

yl3/2 

IF 



= Ko ±— Vff (b,iP) (21) 



where 



1 c 7 ackV 2 m 2 u ml /2 A 2 

K ° ~ 2.947^ e ^ 2 [22) 

16a 196.5 A 2 g 

" 3 mu 24.59 Z3 cm 5 K V2' W 
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C7 = 316.8, a is the Stefan-Boltzmann constant and b = (3/t. The factor of 
2.947 scales the results of Silant'cv & Yakovlev (1980) to agree with the results 
of Cox & Giuli (1968) (for discussion see Hernquist 1985). 

We parameterize the effects of the magnetic field by the anisotropy factor 
for free- free absorption (t]ff). Absorption dominates the opacity through the 
nondegenerate portion of the envelope (Pavlov & Yakovlev 1977, Silant'ev & 
Yakovlev 1980). We use the analytic results of Pavlov & Yakovlev (1977) to 
extrapolate beyond the tabulated values in Silant'ev & Yakovlev (1980); i.e. for 
b > 1000. For b < 1000 we use the results of Silant'ev & Yakovlev (1980). 



4.1 Non-degenerate Structure Equation 

We follow the method of Hernquist & Applegate (1984), but we do not require 
that the conductivity be a power law in p and T. 

In the non-degenerate regime, the thermal structure equation of the envelope 

Z = (24) 
dP g sP K 1 ' 

and we consider an unmagnetized conductivity which is a power law as Equa- 
tion 21 

K = K o — - (25) 

Even in an intense magnetic field, in the nondegenerate regime, the pressure is 
given by the ideal gas law (Blandford & Hernquist 1982) 

P = —pkT. (26) 
m u 

We combine Equations 24-26 with Equation 21, yielding 
dT F 1 fm^V- 1 P"- 1 



dP g s K \Y e kJ Vff (b,ij)T a + s - 1 ' 

As for the structure equations in the degenerate limit, this equation is separable, 
yielding p(T). Because rjff depends on T through the argument b, the relation 
between T and P need not be a power law as in the unmagnetized case. In the 
limit that t]ff = 1, the result of Hernquist & Applegate (1984) obtains. 

More generally, if we take, rjff to be a power law oc b~ 2 (e.g. Tsuruta 
et al. 1972) which is approximately true for b — > oo, we can immediately use 
the results of Hernquist & Applegate (1984) to obtain that the conductivity is 
constant along a solution through the nondegenerate envelope: 

a + S-2 F Y e k 

k = (28) 

a g s m u 

= 7.07 xlO 13 ^ ^ 6 CTg (29) 
-456 9s,u K cm s 
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for a = 2 and 5 = 6.5 as in Equation 21. If we equate this result with our 
assumed conductivity and solve for T, we find that the solution follows 

1/(5-2) / „ \4/9 

(30) 



(\ L/(0-Zj / \ ' 

— ] p «/( 5 - 2 ) = io 6 k(-^) 



where 



2 ( 0m^ y 1 

3 V k J ln6 Ty picai 

PT6 = 71.3 g cm- 3 -=^^A% 2 Z^T;i fi g\% (32) 

V m ^Typical 

and ^Typical is a typical value of (3/t in the envelope, &T yP icai ~ 6 x 10 3 /3. One 
should note that for free- free scattering in the weak- field limit, T oc p 4 / 13 . 

With Equation 30, we can calculate the density at the onset of degeneracy. 
We will assume that at the onset of degeneracy the electron density is given by 
fully degenerate expression and that t«(-1. This yields 

Pnd/d = 3.92 x 10 5 g cm- 3 /?(ln & TyP icai) 1/7 A% 7 Z^T^g-] 17 . (33) 

In principle, electron scattering could also play a role in the nondegenerate 
regime. To simplify the calculation, we neglect its contribution and verify that 
it is indeed negligible. Using Equation 30, we find that the ratio of free-free to 
electron scattering opacity along a solution is given by 

~(F) / \ 14/9 

^ - 1 74 X 10 6 /T 5/9 [ t ) A 4/3 £- 10 /9 T -28/9 7/9 f „ 4) 

~(T) — -1-' 4X 1U Pa R — r "^56 ^26 J eff,6 #s,14> V*^) 

4=0 \V ln ^Typical/ 

where the electron scattering opacity Kp2 1S given by Silant'ev & Yakovlev (1980) . 
Since this ratio increases with decreasing density, we only need to estimate it at 
the maximum density for the solution, i.e. the density at the onset of degeneracy 



~(F) 
K /3=0 



~(T) 
K /3=0 



= 1.36 x 10 3 /3(ln6 Typical )- 6 / 7 4/ 7 z£' 7 T^' 7 g% (35) 

ND/D 



For B — 10 14 G this ratio is greater than one for T e s < 5.9 x 10 6 K, which is 
larger than the effective temperatures considered here. Furthermore, this is a 
conservative estimate of this ratio because generally we cut off the nondegen- 
erate solution where degenerate electrons begin to dominate the heat conduc- 
tion. In unmagnetized envelopes this occurs where the gas is mildly degenerate 
(Hcrnquist & Applegate 1984). We find that this is also the case for strongly 
magnetized envelopes. 
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5 Calculations 



5.1 Strategy 

We have found that the heat transfer equation is not solvable analytically, but it 
is separable in several cases. For the liquid and degenerate region of the envelope 
the solution may be calculated once for each field strength and geometry, and 
scaled to reflect the magnitude of the heat flux and shifted to fit the temperature 
at the low density edge of the region. We can apply this same strategy in the 
solid state only for the case of a purely radial or azimuthal field. Otherwise, for 
the solid region, the temperature as a function of chemical potential will depend 
on the flux and boundary conditions in a more complicated way; consequently, 
the solution must be recalculated for each value of the flux. To match the solu- 
tions across the liquid-solid phase transition and the degenerate-nondegenerate 
interface, we follow the approach of Hernquist & Applegate (1984). At the outer 
boundary, we use the radiative zero solution. 

5.2 Results for the separable structure equations 

As described earlier, in the low-temperature limit when only one Landau level 
is filled the structure equation is simple and may be integrated for a given field 
strength and geometry and the boundary conditions and the dependence on the 
flux may be satisfied after the numerical solution is obtained. 

Given these numerical results, it is straightforward to calculate the core 
temperature for a given surface temperature and flux. However, before fixing the 
boundary conditions, we can note several general features of the results. First, 
for transport along the magnetic field, the envelope becomes nearly isothermal 
at ( — 1 ~ 0.1 regardless of the magnetic field strength. However, for transport 
perpendicular to the field, the temperature rises steadily throughout the range 
of applicability of this formalism. 

Parallel Transport In Figure 1 we present results for the degenerate and 
non-degenerate regime for several magnetic fields with an effective surface tem- 
perature of 10 6 K. In the nondegenerate regime the temperature solution follows 
the power-law given in Equation 30 and the conductivity is nearly constant. In 
the degenerate regime, the conductivity increases dramatically and the temper- 
ature remains nearly constant. For the solution with B = 10 15 G, the discon- 
tinuity in the conductivity at the phase transition is apparent. The results for 
10 14 G qualitatively agree with the results of Hernquist (1985) for this field 
strength. Quantitatively however, we find that the conductivity in the nonde- 
generate regime is thirty percent lower than the earlier result of 10 14 erg/(K cm 
s) and is given by Equation 29. As the magnetic field strength increases, we find 
that the core temperature (or here the temperature at which the first Landau 
level is filled) decreases. This effect results from the increased conductivity in 
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Figure 1: Thermal structure of a strongly magnetized neutron star envelope 
for a radial field. The left panel traces the temperature-density relation with 
B = 10 14 , 10 15 , 10 16 G and effective surface temperature of 10 6 K. The right 
panel traces the conductivity through the envelope. The constant conductivity 
solution appropriate for a purely power-law conductivity law works well through 
the nondegenerate regime. In the left panel, liquid phase exists above the dashed 
curve, and solid phase exists below. 
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the nondegenerate regime where k is approximately proportional to f3 2 and the 
degenerate regime where the quantization of the electron phase space increases 
the conductivity above the zero-field values. 

We take advantage of the simplicity of this analytic technique when calculat- 
ing the thermal structure for hotter and cooler surface temperatures. We do not 
need to reintegrate the structure equations themselves. All that is required is to 
recalculate the boundary conditions at the non-degenerate-degenerate interface 
and the liquid-solid phase transition. Again we find qualitative agreement with 
Hcrnquist (1985). In the nondegenerate regime, the increased or decreased flux 
mimics the effect of changing the field strength depicted in Figure 1 . 

We compare the various results by determining the temperature at the fol- 
lowing densities: p — 1.5 x 10 7 , 4.7 x 10 s and 1.5 x 10 10 g/cm 3 . These are the 
densities at which the lowest Landau level fills for field strengths of 10 14 , 10 15 
and 10 16 G. Moreover, since the matter is nearly isothermal at higher densi- 
ties, these temperatures are close to the core temperature, at least for parallel 
transport. By fitting the results of the calculations, we find that at the lowest 
density: 

T(p = 1.5 x 10 7 g/cm 3 ) oc /T ' 19 ( — J (36) 
and at both the higher densities 

/ f\ 43 

T(p = 4.7 x 10 8 g/cm 3 , 1.5 x 10 10 g/cm 3 ) oc /3" 016 ( - J . (37) 

Figure 2 compares the numerical results with the best-fit power-law relations. 

Perpendicular Transport Modeling the transition between photon and elec- 
tron heat transport is qualitatively different for transport perpendicular to 
the field lines. In the parallel case, the conductivity from electrons typically 
increases rapidly with density, and the transition from photon to electron- 
dominated heat transport is abrupt. For perpendicular transport, the function 
Q decreases with energy, and therefore the conductivity decreases with density 
for fixed temperature. In this case the transition is more subtle. Fortunately, 
the solution does not depend strongly on how this transition is treated, so we 
choose to employ ( — 1 > r to delineate the region where electron conduc- 
tion dominates. The conductivity is not continuous across this transition as is 
apparent in Figure 3. 

We varied the definition of the non-degenerate-degenerate interface and 
found that it had little effect on the T max — T e ff relation. Figure 4 shows how 
the solution changes if we move the interface to a factor of ten higher or lower 
temperature (i.e. ( — 1 > lOr and £ — 1 > r/10). Although near the inter- 
face the solutions differ dramatically, at higher densities the choice has little 
effect. The boundary condition at the transition is unimportant for perpen- 
dicular transport, because the temperature rises quickly with density, and the 
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10" 1G 15 10' 6 -0" 
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Figure 2: The left panel depicts the temperature- flux relation for several field 
strengths and densities. F/g s is given in units of <r(10 6 K) 4 /10 14 cm/s 2 . The 
right panel depicts the tempcrature-magnctic-field relation. The symbols show 
the calculated data points and the lines arc the best-fit power law functions to 
the data. 

solution quickly "forgets" the boundary conditions, in a manner analogous to 
the convergence of the radiative zero solution to the true solution in stellar atmo- 
spheres (e.g. Schwarzschild 1965). This is in contrast to the case where ip ^ tt/2 
where the material quickly becomes isothermal in the degenerate regime. 

We find that for a given effective temperature the core temperature is much 
higher where the heat must travel perpendicularly to the field lines. Further- 
more, we find that for stronger field strengths the effect is more pronounced. 

Effective Temperature Distributions To find the effective temperature as 
a function of angle with respect to the magnetic field we vary the effective tem- 
perature as a function of angle until the temperature where the first Landau 
level fills is constant for the various angles. Unfortunately, where the magnetic 
field is neither radial or tangential we have solutions only in the non-degenerate 
and liquid degenerate regimes. Therefore, for acute angles we must select fluxes 
such that the degenerate solution is entirely in the liquid regime. For the more 
strongly magnetized envelopes we can follow the solution to higher densities; 
consequently, we must use larger effective temperatures for the stronger mag- 
netic fields. Tabic 2 summarizes the parameters for the calculations. Figure 5 
depicts the results for B = 10 14 and 10 16 G. For all but the perpendicular 
case, the envelope has become nearly isothermal by the density where the first 
Landau level fills. 
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T„„=10 5 -° K, g,= 10" cm/s 2 T,,= 10 5 ° K, g,= 1 0" cm/s 2 




2 4 6 8 '0 12 2 4 6 8 10 12 

I°9,q ( P [g/c-n 3 ] ) log. ( p [g/cm 3 ] ) 



T = ,,= 1 5 5 K, g,= 10" cm/s 2 _ . T Iff =10 5 - 5 K, g,= 1C" cm/s 2 

90° lO 1 ^ 00 10 ' 6 




2 4 6 8 10 12 2 4 6 8 10 

log-o ( P [9/ cm3 ] ) lo 9,o ( P [gA^ 3 ] ) 



Figure 3: Same as Figure 1 for the perpendicular case. 



Table 2: Results of the Two Dimensional Calculations 



Field Strength [G] 


T eS {ip = 0) [K] 


T max (C = V20+1)] [Kj 


10 i4 


1.07 x 10 b 


1.12 x 10 8 


10 15 


3.16 x 10 6 


4.79 x 10 8 


10 16 


5.06 x 10 6 


8.13 x 10 8 
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T eff = 1 5 ' 5 K, g s = 1 14 cm/s 2 




2 4 6 8 10 12 

iog 10 ( p [g/cm 3 ] ) 



Figure 4: The dependence of the envelope solution for transport perpendicular 
to the magnetic held upon the dchnition of the non-degenerate-degenerate in- 
terface. We have calculated the location of the interface for (£— l)/r = 0.1, 1, 10 
from left to right. 
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T oIr = 1 6 - K, g s = 1 14 cra/s' T o . f = 1 6 - 7 K, g,= 10 u cm/s 2 




°9,o ( P [9/cm 3 ] ) log 10 ( P [g/cm 3 ] ) 



Figure 5: The left panel shows the temperature structure of the envelope as a 
function of density for B = 10 14 G and T cB = 1.07 x 10 6 K. From top to bottom, 
the results are for tp = 0°, 30°, 60°, 85°, 90° where ijj is the angle between the 
magnetic field and the radial direction. The right panel depicts the temperature 
structure for B — 10 16 G and T e g — 5.06 x 10 6 K. The solutions are constrained 
to have the same temperature at the density where the first Landau level fills 
(denoted by the bold circle). In both panels, liquid phase exists above the 
dashed curve, and solid phase exists below. 
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Figure 6: Flux as a function of angle for B = 10 14 , 10 15 and 10 16 G. The solid 
curve is cos 2 i[>. 
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Figure 6 shows the flux as a function of angle for all of the two-dimensional 
calculations. The agreement between the flux distribution and a simple cos 2 ip 
law is striking. Greenstein & Hartke (1983) have argued that if the conductivity 
is constant through the envelope, the flux will follow a distribution of the form 
A cos 2 ip + B sin 2 tp. Although in the nondegenerate regime the conductivity 
along a T(p) solution is nearly constant, in the degenerate regime it varies by 
several orders of magnitude. Furthermore, the nondegenerate layers do not 
throttle the heat flux; if they did, one would expect little variability, as the 
conductivity parallel and perpendicular to the field are nearly equal in the strong 
field limit. 

We look to the degenerate structure equation for the liquid state to explain 
the remarkable agreement with a cos 2 ip distribution. ^From examination of 
Equation 18, we see that if the conductivity transverse to the field is neglected, 
we can make the replacement 

F^Fcos 2 ip (38) 

and recover the thermal structure equation for ip = 0. We determine where this 
approximation is valid by comparing the transverse and parallel components of 
the conductivity tensor 

K yy,ei _ Z 2 a 4 2 QeijC'-, P) _o W 1n -ll -^26 ^ QeijC'jP) /o Q n 

^~4^ 2 /33 (C 1} ^(C;/3) ™ 8 - 22xl ° A 2 6 ^ 4 % e4 (C;/3)- m 

At first glance, it appears that the transverse conductivity is negligible through- 
out the degenerate regime. However, the functions Qei and 4> e i complicate the 
discussion. Specifically, 4> ei — > and Q ei — > oo as C — » 1; therefore, transverse 
conduction is likely to be important in the nonrelativistic portion of the degen- 
erate envelope. Figure 5 shows that this is indeed the case. For B = 10 14 G, the 
solutions for ip < tt/2 are nearly identical for p > 10 6 5 g/cm 3 or £ > 1.1. As 
£ approaches unity, the ratio of the conductivities increases without bound, the 
transverse conductivity may no longer be neglected, and the runs of temperature 
with density begin to diverge. 

Empirically, we find that in the region of the envelope which most effectively 
throttles the flux, the transverse conductivity may be neglected for ip < tt/2 
without introducing significant error. 

5.3 Observed Flux Distribution 

We follow the technique outlined by Page (1995) to calculate the observed fluxes. 
However, unlike Page (1995) we evaluate the double integrals over the visible 
portion of the neutron star surface directly. We use the cos 2 ip rule to calculate 
the photon distribution function at the surface, and so we do not define a grid 
of precalculated distribution functions as Page (1995) does. 
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As a first approximation, we focus on the variation of the observed bolomctric 
flux with the angle p of the line of sight with the magnetic dipolc axis. This 
angle is a function of the phase angle (7), the inclination of the dipole to the 
rotation axis (a) and the line of sight to the rotation axis (() 

cos p = cos £ cos a + sin £ sin a cos 7 (40) 

(Greenstein & Hartke 1983). For simplicity, in the discussion that follows we 
will take a = ( = tt/2; therefore p = 7 and we refer to ip = 0,7r as on-phasc 
and ip = it/2, 3tt/2 as off-phase. 

We repeat the calculation for several values of the stellar radius (with fixed 
mass) to determine the effects of general relativity on the light curves: gravita- 
tional redshift and the deflection of null geodesies (self-lensing, or more concisely 
"lensing"). 

To quantify the effect of gravitational lensing on the light curves of magne- 
tized neutron stars, we calculate the mean value of the bolometric flux emitted 
by the surface over the visible region of the star. We assume that the flux at 
a given location on the surface is proportional to cos 2 tp where ip is the angle 
between the radial direction and the magnetic field. 

For clarity, we treat the gravitational redshift separately. Figure 7 depicts 
the ratio of the mean value of the flux over the visible portion of the star to 
the flux that would be emitted if the magnetic field were normal to the surface 
throughout (i.e. an isotropic temperature distribution). 

In the limit of infinite radius, i.e. if lensing is unimportant, we find that for 
a uniform field 

f(p = 0) = \f{ip = 0) and f(p = tt/2) = = 0). (41) 

For a dipole field, the calculation is slightly more complicated. First, we used 
Equation 10 to determine the angle of the field with respect to the radial direc- 
tion. Secondly from equations 36 and 37, we find that the emergent flux is a 
function of the field strength. For a dipole configuration, the magnitude of the 
field varies as 

oc \/3cos 2 (9 + 1 (42) 

along the surface of the star. Since we are most interested in fixing the internal 
temperature at high densities we assume that the flux is proportional to B 0A 
from Equation 37 which reduces the flux for ~ 7r/2 further beyond the cos 2 ip 
rule. We obtain 

f(p = 0) = 0.663/(V> = 0) and f(p = tt/2) = 0.393/(V> = 0). (43) 

If we did not include the effect that the flux is a function of field strength as well 
as orientation we would have obtained 0.717 and 0.444 for the above values. 

We find that the mean flux is greater for the dipole configuration than for 
a uniform field for all viewing angles if R/R s < 5, and that the variation in 
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Figure 7: The average of the bolometric flux over the visible portion of the 
neutron star for ip = 0, ir/2. The upper pair is for a dipole, and the lower is a 
uniform field configuration. 
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the light curve is generally smaller. We have taken Mns = 1.4M Q , yielding a 
Schwarzschild radius, R s , of 4.125 km. The theoretical predictions for the radius 
of a 1.4M© neutron star range from 6.5 km to 14 km (Thorsson, Prakash & Lat- 
timer 1994; Pandharipande 1971; Wiringa, Fiks & Fabrocini 1988), depending 
on the details of the equation of star at supernuclear densities. 

As Page (1995) found, lensing dramatically reduces the variation of the ob- 
served flux with phase by making more than half the surface visible at any 
time. Interestingly, for the range of radii 7.248 - 8.853 km, we find that the 
flux is greater when the magnetic poles are located perpendicular to the line of 
sight. For this range of radii over 90% of the surface is visible. Emission from 
both of the hotspots reaches the observer leading to a larger flux. For radii less 
than 7.248 km the entire surface is visible and again the peaks are on phase. 
Page (1995) found a similar effect for the same range of radii. 

The emitted spectra from the visible portion of the neutron-star surface is 
the sum of blackbody spectra of various temperatures. To determine the emitted 
spectra, we calculate the distribution of blackbody temperatures on the surface; 
i.e. we estimate the distribution function d//dT e ff. With this distribution func- 
tion, it is straightforward to calculate the emergent spectrum averaged over the 
visible portion of the surface as 

U ~J dT cS aT c %A^eMrito/kT)-l dIcS - {U > 

The factor of (tT* s converts the flux to an effective area of emission (A). The 
calculation of d//dT e fj is numerically more tractable than dA/dT e g and allows 
us to account for the total energy emitted more reliably. 

We calculate df/dT c g in similar fashion to /. To expedite the calculation, 
we note that given the cos 2 ip rule the neutron star surface has a limited range 
of effective temperatures, specifically between and T e s(ip — 0). Consequently 
we define 

For observations on-phase (ip = 0) and without lensing (R — > oo), the flux- 
weighted temperature distribution can be calculated directly if the envelope is 
uniformly magnetized. It is given by 

^=4T 7 /(V = 0). (47) 

The result for the dipole cannot be written explicitly and is not illustrative. 
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For a general geometry (</? ^ 0), we expand this function in an orthonormal 
basis on the interval to 1. Specifically, we assume that for < T < 1 

^ = f>Qi(T) (48) 
ai 1=0 

and zero otherwise, where 



Ql (T) = V2lTTPi(2T - 1). (49) 

The Pi(x) are the Legendre polynomials. From the properties of these orthonor- 
mal functions, we have A = /, and Ai is calculated by inserting the weighting 
function Qi(T) into the integrands in the calculation for /. We recall that we 
have assumed, T — cos 1 / 2 ip 

Using an orthonormal basis dramatically speeds the calculation of the dis- 
tribution. Additionally, because the Pi{x) are polynomials, it is straightforward 
to calculate conventional statistics of the distribution 

= if 1 f^df^^ + 1 - (50) 

/ = 30^ + — ^ + 3' (51) 

Unfortunately, with this basis is impossible to insist that distribution is every- 
where non-negative, i.e. that no temperatures contribute negative flux. How- 
ever, if a sufficiently large number of Ai are calculated, the intervals where 
d//dT < can be made to be arbitrarily small and to have an arbitrarily 
small contribution to the total flux. We compare the results of the expansion 
with Equation 47 and find that the maximum relative error in the expansion 
coefficients between the two techniques is approximately 9 x 10~ 5 . 

Figure 8 depicts the results of this calculation for four values of the stellar 
radius with ip — 0,7r/2. In the left panel, we see in the absence of general rel- 
ativistic effects that when the neutron star is off-phase more flux is produced 
at lower blackbody temperatures than when the magnetic dipole is pointing to- 
ward the observer. For the smallest radius considered (R — 6.25 km), the entire 
surface of the neutron star is visible and a large portion of the front hemisphere 
has a second image. In this case, both the flux-weighted temperature distri- 
butions at ip — and tt/2 are peaked at the maximum effective temperature. 
However, the distribution off-phase has a more well populated tail extending 
toward lower temperatures than on-phase. 

The right panel depicts a value of the stellar radius (R — 8.85 km) where 
there is practically no variation of / with phase. Additionally, we see that the 
flux-weighted temperature distributions are nearly constant with phase. Also 
depicted is the temperature distributions for R = 7.9 km, the radius where the 
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Figure 8: The fractional distribution of observed flux as function of the surface 
blackbody temperature. The left panel depicts the distribution for the minimal 
(R/R s = 1.52) and maximal radii (R = oo) considered. The right panel depicts 
the distribution at a radii where f((p = 0) = f(cp = tt/2) (R = 8.85 km) and 
where the off-phase peaks are maximized (R = 7.9 km). All the curves are 
normalized to have an integral of unity from T = to 1 



off-phase peaks are maximized. Here again, the distributions do not change 
appreciably with phase. We conclude that for 7 <J R <, 9 km, it would be 
difficult to detect variation of the spectra with phase, if the neutron star surface 
indeed radiates as a blackbody. 

To calculate the spectra themselves, it is convenient to define the following 
functions, 

15 „„ f 1 „ 1 1 



?m = ^&l qi(t)^ ; dT, (52) 

n Jo J- exp(w/J J — 1 

where u> — Tvuj / kT e s(ij) — Q), so that 

Since Ai has units of flux, we obtain the correct units of flux-time for / w . 

To convert to the observed spectra, we must account for gravitational red- 
shift and interstellar absorption; we obtain (Page 1995) 

/° bscrvcd M - jl(e- A ^) e - A = *~ e -"H*(«) (54) 

where 

R^ ee Re As (55) 
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where e s is given in Equation 6. The final term accounts for interstellar 
absorption and D is the distance to the neutron star. For <r(E) we use the 
Morrison & McCammon (1983) cross-sections. 

Figure 9 depicts the spectra in the ROSAT energy range for two neutron 
stars. Each of the spectra (light curves) is well fitted by a blackbody (heavy 
curves) at the mean effective temperature with an additional hard component. 
The mean effective temperature (T mean ) is defined by 



for Nh = 0, i.e. it is the equivalent blackbody temperature that accounts for 
all of the energy emitted from the neutron star surface. The hard component is 
most significant when the star is observed at right angles to the magnetic axis, 
and originates from the portions of the hot polar caps that are visible even when 
the star is off-phase. 

In Appendix A, we present two XSPEC models which are available over the 
WWW. With these models, one may simulate observations from various x-ray 
instruments to estimate the observed pulsed fractions for the models discussed 
in this section. We give an example in Figure 11. 

5.4 Neutron Star Cooling 

We can use the results of the preceding sections to determine the effects of 
the magnetic field on neutron star cooling rates. Specifically, we take the ra- 
tio of the total flux from the surface with and without a magnetic field for 
the same core temperature. We have used the results of Hernquist & Apple- 
gate (1984). To determine the core temperature for a given flux we combine 
their equations (4.7) and (4.8), switching from the first relation to the second 
when the surface effective temperature drops below 4.25 x 10 5 K. The results 
do not depend qualitatively on whether equation (4.7) or (4.8) of Hernquist & 
Applegate (1984) is used. 

To compare our calculated temperature-flux relation with the results for 
isotropic heat transport, we multiply the fluxes for the magnetized envelopes 
by 0.4765 to account for a dipole field configuration. Additionally, we assume 
in our calculations that the envelope is isothermal above the density (p max ) at 
which the first Landau level fills and use the temperature (T max ) at this density 
to estimate the flux in the unmagnetized case. Only for the strongest field 
strength considered (B — 10 16 G) do our analytic calculations extend to the 
core density assumed by Hernquist & Applegate (1984) of 10 10 g/cm 3 ; however, 
we do not expect the results to be strongly sensitive to this cutoff density as 
our solutions (for ip ^ 7r/2) are nearly isothermal at high densities. 

Figure 10 depicts the results of this comparison. We find that for the weakest 
field strength considered (B = 10 14 G), the magnetic field has little effect on the 
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Figure 9: The observed spectra from a neutron star with T e g(ip = 0) = 7.5 x 
10 5 K from a distance of 250 pc with an intervening absorption column of 
Nh — 10 20 cm -2 . The light curves show the spectra, and the heavy curves show 
blackbody spectra at the mean effective temperature of the neutron star. For 
the R = oo model, we have taken R/R s — > oo while R — 20 km to give the 
surface area of the neutron star, while neglecting general relativistic effects. 
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Figure 10: The ratio of fluxes in magnetized envelopes to those in unmagnetized 
ones. We have assumed that the envelope is isothermal above the densities given 
in the legend to estimate the unmagnetized fluxes. 
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total luminosity of the star. However, for cooler core temperatures and stronger 
magnetic fields, the difference in the luminosities can be up to a factor of ten. 
The flux ratio is sensitive to the core temperature because Equations 36 and 37 
have a slightly different power-law index than the model assumed for the un- 
magnetized envelope (0.392). The inflection in each of the curves occurs when 
the material near non-degenerate-degenerate interface melts as the core temper- 
ature increases. All the curves swing upward for high values of T max because for 
high fluxes (i.e. high core temperatures) our assumption that the temperature 
is constant for p > p max no longer holds. 

Because strongly magnetized neutron stars emit significantly more flux, we 
expect that the thermal history of magnetars should be dramatically different 
from that of neutron stars with weaker magnetic fields. We discuss this issue 
further in Heyl & Hernquist (1997). 

6 Conclusions 

We have presented an analytic technique for calculating the thermal structure 
of ultramagnetized neutron star envelopes. We use the exact thermal conduc- 
tivities in an intense magnetic field of Silant'ev & Yakovlev (1980) and Hern- 
quist (1984) in the non-degenerate and degenerate regimes, respectively. We 
make two simplifying approximations. We assume that the interface between 
degenerate and non-degenerate material is abrupt. Hernquist (1985) numeri- 
cally calculated the thermal structure for B = 10 14 G without this assumption. 
Our agreement with this earlier result shows that an abrupt interface is a good 
approximation. Secondly, we use a standard simplification in the study of stellar 
atmospheres which is to use the radiative zero solution to fix the outer boundary 
condition (Schwarzschild 1965). Because the equation for the thermal structure 
in the outermost layers is qualitatively similar to the relation for an unmag- 
netized envelope, we conclude that as in the unmagnetized case (Hernquist & 
Applegate 1984), this is also an accurate approximation. 

A distinct approach is to treat the entire problem numerically, which allows 
us to estimate directly the possible errors that our simplifying assumptions 
introduce and alleviates the problem of how to treat the regions where more 
than one Landau level is filled. Numerical models for parallel and transverse 
conduction for B = 10 12 — 10 14 G and p < 10 10 g/cm~ 3 are presented in Heyl 
& Hernquist (1998). 

We find that the relation between transmitted flux, core temperature and 
field strength may be approximated by a power law and that the effective tem- 
perature is proportional to cos 1 / 2 ip where tp is the angle between the radial 
direction and the local direction of the magnetic field. Using the geometric re- 
sult, we calculate the observed spectra as a function of viewing angle including 
the effects of general relativity for dipole and uniform field configurations. We 
extend the conclusions of previous work. If the surface is assumed to radiate as a 
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blackbody and neutron stars have radii within the currently accepted range, the 
anisotropic heat transport induced by a dipole field configuration is insufficient 
to produce the observed pulsed fractions even for ultramagnctizcd envelopes. 

Pavlov et al. (1994) argue that in addition to the transmission of heat 
through the envelope, the emission at the surface is also strongly anisotropic. 
Anisotropic emission can naturally produce large pulsed fractions even when the 
temperature on the surface is uniform (Zavlin ct al. 1995, Shibanov et al. 1995). 
Additionally, the composition of the atmosphere may have a profound effect on 
the emergent radiation. A magnetized iron atmosphere produces substantial 
limb darkening for ip ~ 7r/2, and the decrement is strongest at high energies 
(Rajagopal, Romani & Miller 1997). It is straightforward to graft these atmo- 
spheres onto the thermal envelopes calculated here to obtained the observed 
time-dependent spectra for a variety of realistic neutron star models. These 
effects along with anisotropic conduction may be sufficient to account for the 
large observed pulsed fractions. 
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Rather than present results for specific instruments and band passes, we supply 
our results in machine-readable form. We have calculated neutron star spectra 
for several values of R/R s , and <p = 0, n/2. We assume a dipole field configura- 
tion and the cos 2 tp rule. The model is calculated for 



where R s = 2GM/c 2 . T c s t00 (tp = 0) may be varied by applying a redshift 
and renormalization to the spectra. For ease of use by the x-ray astronomy 
community, we have created XSPEC table models with the data. 

Within XSPEC, the models may be convolved with the response matrix for 
various x-ray instruments and compared with observed spectra. By using the 
redshift (z) and normalization (K) parameters we may obtain models for dif- 
ferent effective temperatures and radii: 
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where i?oo,km is the source radius in km as observed at infinity, i.e. 

R 00 =Re A «=R\l-^-\ (60) 

from Equation 55 and D\q is the distance to the neutron star in units of 10 kpc. 
This particular choice is consistent with the definition of the bbodyrad model 
in XSPEC. 

Each additive model contains the single interpolation parameter R s /R which 
ranges from to 0.6601. As an illustration Figure 11 depicts one of the models 
convolved with the ROSAT PSPC response function for R/R s — > oo, R = 20 
km and D — 250 pc with Nh — 10 20 cm~ 2 . Here, we have used the wabs 
model to calculate the interstellar absorption which assumes the Morrison & 
McCammon (1983) cross-sections. 

The errorbars are calculated for an exposure of 10 4 seconds. For these pa- 
rameters, the variation of the thermal flux with phase is apparent in the spectra. 
However, as we saw in the previous sections as R/R s decreases, the variation in 
the thermal flux weakens. 

The table models are available at the following URLs: 
http : //www, ceo . caltech. edu/^jsheyl/analytic_ns/pO .fit 
http : //www . ceo . caltech. edu/^j sheyl/analytic_ns/p90 .fit 

The XSPEC software itself is available at 

ftp : / /legacy . gsf c . nasa. gov/ sof tware/xanadu/ , 
and an online manual is provided at 
http : //www .merate .mi . astro . it/~xanadu/xspec/u_mamial .html 



for ip — 
for (p = 7r/2. 
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Figure 11: The calculated model spectra convolved with the ROSAT PSPC 
response matrix for ip = (upper points) and tp = n/2 (lower points). The 
model parameters are described in the text. The errorbars are for an exposure 
of 10 4 s. 
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